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The details of fluctuation free energy conservation in the gyrofluid model are examined. The 
. . . polarisation equation relates ExB flow and eddy energy to combinations of the potential and the 

' density and perpendicular temperature. These determine the combinations which must appear 

' under derivatives in the moment equations so that not only thermal free energy but its combination 

, with the ExB energy is properly conserved by parallel and perpendicular compressional effects. The 

resulting system exhibits the same qualitative energy transfer properties as corresponding Braginskii 
' or Landau fluid models. One clear result is that the numerical model built on these equations 

, is well behaved for arbitrarily large perpendicular wavenumber, allowing exploration of two scale 

phenomena linking dynamics at the ion and electron gyroradii. When the numerical formulation is 
^■f^ ' done in the globally consistent flux tube model, the results with adiabatic electrons are consistent 

, with the "Cyclone Base Case" results of gyrokinetic models. 
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PACS numbers: 52.65.Tt, 52.35.Ra, 52.30.-q, 52.25.Fi 



^ ■ I. INTRODUCTION - GYROFLUID ENERGY CONSERVATION IN GENERAL 

Gyrofluid modelSjWhose most prominent application has been to tokamak core turbulence as exemplified by the 
■ Cyclone Base Case were originally constructed to incorporate finite ion gyroradius effects at arbitrary order into 
, ^ ' simple computations of turbulence occurring in largely two dimensional fluid experiments [3] . To treat ion temperature 
^ • gradient (ITG) turbulence the temperatures were incorporated and the model acquired several new advection terms, 
producing nonlinearities as well as drift frequency corrections, resulting from the effect of temperature fluctuations 
on the gyroaveraging operator However, although the nonlinearities were incorporated, much of the analysis 

involved linear frequencies and growth rates and the nonlinearities were added largely as an afterthought. Specifically, 
there was no complete energetic analysis of the type familiar from drift wave turbulence work 0, 0, 13 • The GEM 
model was introduced previously in the context of energetic considerations, including the correspondence between fluid 
drift and gyrofluid models under drift ordering 0|. We develop the energetics for the six- moment model including 
temperature and parallel heat flux dynamics herein, placing this model on a secure energetic footing for the first time. 
^-H Under drift ordering |10| , this energetics involves "fluctuation free energy," in which the thermal free energy enters as 

T— I the average squared amplitude of the density fluctuations Q , with additional contributions from the average squared 
amplitude of the temperature fluctuations in the appropriately generalised models Q- The rest of the free energy is 
made up of contributions due to the ExB energy involving the electrostatic potential, the magnetic energy involving 
the parallel magnetic potential, and the parallel free energy in not only the parallel velocities but also the parallel 



heat fluxes . A properly constructed model should conserve this free energy in all processes except those involving 
, ^ clearly identifiable sources (gradients) and sinks (dissipative processes such as resistivity, thermal conduction, or 
^ ' Landau damping). Many models neglect this consideration because small errors in the energetics lead simply to 



negligible contributions to the growth rate in a linear model. However, if the model is to be useful in a turbulent 
setting, the energetics must be consistent in order to achieve a reliable saturated state in which the salient energetic 
processes and the turbulent transport can be statistically measured. 
. ^ ' T^hB processes linking the thermal free energy to the ExB turbulence are of physical importance because the free 
k> I energy sources and sinks are not in the perpendicular equation of motion. The energy source given by the general profile 
5—1 ' gradient, is in the equations for the thermal state variables (density, temperatures). The dissipation processes are in 
^ the equations for the parallel flux variables (current, heat fluxes). In saturation, the ExB energy itself is maintained 
as a statistical balance between various conservative transfer effects which connect to these sources and sinks in the 
other parts of the dynamics (this neglects certain rotation damping processes which are often not considered). In this 
balance it is important that the conservative nature of such processes such as shear Alfven dynamics or interchange 
effects is maintained in the model. The simplest example is an isothermal two dimensional magnetohydrodynamical 
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Vi0 = -TMne) (1) 



(MHD) interchange model, given by 

dfi — 

— f + VB • V (ne + ne) = ne/C(0) (2) 
ot 

where the perpendicular Laplacian and curvature operator are defined by 

vi -V • [bx(bxV)] (3) 



' (BxV) (4) 



-B2 

respectively, with B — Bh the equilibrium magnetic field. The ExB velocity is given by 

V£ = 4bxV0 (5) 

The curvature terms in Eqs. H1I2|) respectively represent quasistatic compression of the diamagnetic current and the 
ExB velocity. 

Under drift ordering 10], the background parameters are constants except where operated upon by v^; • V, the 
ExB advection. The magnetic field is treated as constant except for the existence of /C(). In the advection term, the 
velocity is treated as divergence free; the finite ExB divergence is accounted for by JC{(j)). If we multiply these 
equations by — </> and {Te/ne)ne, respectively, and integrate over the entire spatial domain, we find, neglecting surface 
terms. 



dt 2 52 



Vi0 



= — / dATeUelC 



(0) (6) 



^y^A^ (^^^ ^-J rfATeneVB • Vl0gne+ JdATeUeJC (7) 

where / dA by itself gives the total volume. These two lines give the evolution of the ExB drift and thermal free energy, 
respectively. The interchange effect represented by /C() transfers free energy between these two pieces conservatively, 
as we would expect from a compressional process, because /C() is at once a total divergence and a first order linear 
differential operator. The source is given by the advection of the background gradient, proportional to the average 
flux. This model will saturate only if there is a loss process at the boundary, or some nonlinear effect enters to cause 
the source to go to zero, or some additional effect is considered in the model by which an explicit sink term appears 
to balance the source (a detailed analysis of how this functions in a three dimensional drift Alfven model is given in 
Ref. i). 

In a gyrofluid model, there is no equation for the vorticity explicitly involving time derivatives. Instead, at the same 
level of sophistication as above, there is an evolution equation for each of the species' gyrocenter densities, as opposed 
to space densities. In each density equation, the potential is gyroaveraged using a suitable convolution operator which 
is described by a kernel in Fourier space: 



G 



W ^E^k^'^k,^''^"" (8) 
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The usual form for Gk^ is Fg (bi), which is an average of the single particle form jQ{k±v±/fli) averaged over the 
perturbed distribution function 0]. The argument is bi = k\p1, where pi is the thermal gyroradius and Q,i is the 
gyrofrequency. Unless phenomena below the ion scale pi are considered, electron gyroradius effects are usually ignored, 
so that Gk^ = 1 for them. The potential is given by a polarisation equation which equates the two space densities, each 
given by a combination of the gyrocenter and polarisation densities, the latter involving the gyroscreened potential. 
The equations appear as 

+ vb • V (ue + rie) = rie/C \ 4>j ^/C {rie) (9) 
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^ + ub • V {rii + rii) = niK. (^^g) + (?^^) (10) 

^ + fp?Vi0=^ (11) 

Gyroscreening is distinct from gyroaveraging; in general this is given by another operator whose usual kernel in 
kj^-space is Ti^ibi) — 1. The low-Zc^ form of this is just — 6j, which in real space yields the p^V^ form used in Eg. 1111 
It is important to note however that the same operator is involved in the gyroaveraging of the potential (Eq. (SJ as 
in the conversion of the gyrocenter density, n^, to the ion space density given by the left side of Eq. Ijllfl . For the 
ions, the operator is G; for the electrons, it is unity. Consistent with this, the ExB advection of the ion density occurs 
with the gyroaveraged potential, with velocity 

UB--|bxV(^G (12) 
r> 

while the electrons are advected by the "bare" version, the same v^; as in Eq. (0). The curvature operator acts on 
the total force potential of each species, respectively, i.e., the diamagnetic velocity divergences are necessarily kept in 
the model, which cannot take the MHD form if both pressures are to contribute to charge separation. 

The point to be made here is the way in which the polarisation equation, Eq. (|11|1 . determines the acceptable form 
for many of the differential operators in the moment equations. The ExB energy and the electron and ion thermal 
free energies are given by 



2 B 

respectively. Using Eq. Hll() and the Hermitian property of G, we may recast [/^ as 

UE = e— e— (14) 

which shows the ion and electron contributions separately. By similar means, the time derivatives are given by 
OUe ~ drii ~dn 



dt 



^^'^'W ~ ^'^'W ^ ^ J d-^^enelC - J dATiU.IC (^4>g^ (15) 
dAT.n.v..Vlog.. + /.AT.?i./c(?) (16) 



dUe 

dt 



^ = - J dAT,n,UE • Vlogn, + J dKT.n.JC ((/-g) (17) 

From this we can see that, since both polarisation and thermal energy for each species now follow from its density 
equation, the density and corresponding gyroaveraged potential must appear together under spatial derivatives in all 
conservative processes in order for the energetics to remain consistent. For the ions this is (pc + {Ti/nie)ni, while for 
the electrons it is </) — {Tf,/nf,e)nf,. Thus, for the ions, the combination (pQ + {Ti/nie)ni multiplied by the curvature 
term in Eq. IjlUI) yields a term which vanishes under the integration, conserving the free energy, and the same occurs 
for the electrons upon multiplying the curvature term in Eq. @ by — {Te/nee)ne- 

This is relatively trivial for such a one-moment gyrofluid model, but the central point is clear: the gyrofluid closure 
arising from gyroaveraging must be done the same way in the polarisation and in the moment equations. That is, 
in this case, the operator G in the gyroaveraged potential for a given species must be the same as the one operating 
upon the corresponding density in the polarisation equation. In a model which incorporates the temperatures, this 
extends to another gyroaveraging operator acting upon the temperature, and the same one acting upon the potential 
producing extra finite gyroradius advection terms, plus corrections to the temperature wherever the latter appears 
under parallel gradients or curvature terms. This is much more involved and will be treated in the next two sections, 
one discussing the problems with the currently standard gyrofluid model, and the next one formulating the GEM 
model. Then, the last section shows the applications of GEM to the damping of kinetic shear Alfven waves, to the 
hyperfine electron gyroradius scale turbulence problem, and to the Cyclone Base Case which provides the standard 
benchmark. 
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II. DRIFT ORDERING AND NORMALISATION CONVENTIONS 



With many constant factors containing the physical units to be carried about while manipulating the equations, 
it is convenient to go to a system of normalised units. While this is somewhat arbitrary, two examples are useful in 
elucidating the salient units for gyrofluid turbulence. One of these is the prospect of force balance in the electrons 
along the magnetic field. There are several effects which affect the evolution of the electric current in the electron 
Ohm's law, but these all act to mediate the response of the electrons to the two static parallel forces: the pressure 
gradient and the static part of the parallel electric field, given by 

/eV||0 \ , ^ 

rieCV II - V||Pe = Pe (^7 V||logPel (18) 

assuming small disturbances from the equilibrium, and a field line geometry in which the parallel gradient for finite 
sized disturbances is not allowed to vanish (see below), a quasistatic balance between these two forces yields the 
following relationship between the disturbances: 

e Pe 

If the parallel responses also equalises the temperature along the field lines, we then have the combination 

1 = (20) 

This situation is called adiabatic electrons, and the response to the force in Eq. H18|l is called the adiabatic response. 
The importance of the pressure force is what departs gradient driven turbulence in general from the world of MHD. 
The adiabatic response indicates the useful units for the electrostatic potential and all the thermodynamic state 
variables, which will be scaled in terms of Te/e, rie, or Tg following the forms in Eq. (|20|l . 

The other useful example is closer to the idea of the gyrofluid formulation: the relationship between ion inertia and 
polarisation. In the low-fcj^ limit, the ExB energy is the same as the fluid one, 



Ue = n,M, ^ - ' ' 



(21) 



2 2 B 

given in Eqs. H6I13|I . The functional derivative Ue with respect to (p leads to the polarisation density 

given in Eq. Hll|) . Expressing (j) in terms of Tg/e and the densities in terms of rige, we find 

= —pI^It^ (23) 

where ps is the drift scale given by 



e2S2 

This is equally well described in terms of the ion gyroradius, 

— = -7^P?Vi-^ (25) 

The role of the drift scale remains, however, even if Ti — > 0, because it is a purely inertial phenomenon. With several 
species present, Te and ps make good choices for the basis of the normalisation scheme, and we will use them herein. 

We work in terms of an arbitrary number of charged fluids, each with a particular background density and temper- 
ature, and mass and charge state per particle. In terms of the electron density and temperature, rig and Tei the mass 
of a main ion species Af^, and the unit of charge, e, we have a normalised background charge density 

ttz = (26) 
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temperature to charge ratio 



and mass to charge ratio 



T 

(27) 



ZT, 



/i, ^ (28) 
^ ZMi ^ ' 

for each species labelled by z. For electrons, these are Og = Tg = —1 and = —me/Mi. For each species, the 
background mass density is given by Uzfiz, pressure by QzTz, and squared gyroradius by pi = /JzTz, in units of 
UeMi, Pf, — rieTe, and p^, respectively. In fusion applications the main ion is often considered to be deuterium, 
and so the latter mass Mo is used for Mi. However, with discharge experiments run at fusion-relevant normalised 
parameters using several disparate ion types becoming increasingly important it is important not to make this 
choice automatic. 

We work under drift ordering, also called gyrokinetic ordering. There are two basic perpendicular length scales: 
the drift scale ps and the background profile scale L±. Drift ordering assumes their ratio to be small, 

6 - ps/L^ « 1 (29) 

where S is called the drift parameter. It also assumes that the relative amplitude of all the disturbances is small by 
this same order, for example, 

and similarly for the parallel flux variables (velocities, heat fluxes) in terms of the sound speed Cs given by 

= Te/Af, (31) 
It also takes a "maximal" ordering with respect to perpendicular wavenumbers and the drift scale, 

k^Ps - 1 (32) 
while assuming a flute/drift ordering for the parallel wavcnumber, 

fc||/fci=0(J) (33) 

The significance of these statements taken together is that the dynamics can be very nonlinear, in the sense that 

• V '--^ d/dt, even though the disturbance amplitude remains small. It therefore follows that all nonlinearities are 
dropped except for the quadratic ones represented by ExB advection (v^; • V in a fluid model, additionally all its 
finite gyroradius relatives in a gyrofluid model) and, in an electromagnetic model, the "magnetic flutter" nonlinearities 
represented by the contributions due to the magnetic disturbances in the parallel gradient V||. With the flute ordering 
we assume that the disturbed and undisturbed parallel gradient pieces are of similar size. 

An important implication of drift ordering on the treatment of the geometry concerns the divergences of the various 
perpendicular drift fluxes. The ExB velocity in an inhomogeneous magnetic field has a finite divergence, so that both 

• Vn and nV • would enter a fiuid model. However, under drift ordering, both the profile and disturbance are 
advected at the same order by v^;, while the factor of n multiplying the divergence is treated as a constant parameter. 
This forces an explicit split between the advection and divergence terms. Due to their various cancellations 0,^3, 
the diamagnetic fluxes enter only in the divergences, so this splitting concerns solely the ExB velocity in a fluid 
model, plus the associated finite gyroradius forms in a gyrofiuid model. Under drift ordering, then, the advection 
term appears as a pure Poisson bracket form (between the two perpendicular coordinates in a field aligned treatment, 
or between all three pairs in general), multiplied by a constant coeflticient, and the divergence term appears as a 
curvature term (through /C, as above), whose properties are that (1) it is a pure divergence, (2) it is a first order 
differential operator with the linearity property, and (3) all curvature terms are linear in the dependent variables. 

The further impact of drift ordering on the treatment of the magnetic geometry is summarised in the statement 
that if the three coordinates are aligned to the magnetic field such that one of them is parallel (s), one is radial {x, 
across flux surfaces, down the gradient), and the third (y) has vanishing projections to both the equilibrium magnetic 
field and the background gradient, then the variation of the geometry is retained only in s. Field aligning means 
that only one contravariant component of B is nonvanishing, in this case B'^ . These are the basic statements of field 
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aligned coordinates and magnetic flux tube geometry, as explained elsewhere 0, 0| . The principal consequence is 
that the perpendicular Laplacian involves only x and y, and the metric coefficients depend only upon s. Note 
especially that this includes the variation of the field strength B, whose variation along the magnetic field is retained 
but commutes with V^^ (hence, appears as a coefficient in the normalised gyroradius in all gyroaveraging and 
gyroscreening operations). Although the equations of the standard gyrofluid model and of GEM are expressed in 
covariant terms, the flux tube geometry enters when they are discretised in a numerical representation, and their 
essential coordinate dependence is reflected in the abovementioned split between advection by and divergence of the 
drift velocity. We will save further comment on this for the sections (below) describing the numerical scheme. 

Thermodynamic state variables are normalised in terms of their background quantities, the electrostatic potential 
in terms of Te/e, and the parallel magnetic potential in terms of i?oPs- Additionally, a factor of 5 is folded into all 
the normalisations, leaving the ExB advective nonlinearities coefficientless. This is expressed as 



for the potentials, and 



Te BoPsPe 



5-'- Tz^S-'^ (35) 



n^TzCs 



(36) 



for the state and parallel flux variables respectively, where 



is the electron dynamical beta (this enters through Ampere's law). The density and temperature are therefore treated 
on an equal footing with respect to their flux variables, the velocity and heat flux, respectively. Additionally, in 
both GEM and the standard gyrofluid model, parallel and perpendicular temperatures and parallel-parallel and perp- 
parallel heat fluxes are treated separately, leaving a six moment model, for both ions and electrons as previously [Iq , 
and for all species herein. 

We leave Vy normalised in terms of L±. Hence, the size of the contravariant magnetic unit vector component f 
is comparable to L±/qR, where q is the magnetic pitch parameter (toroidal/poloidal magnetic fleld, contravariant 
component ratio), R is the toroidal major radius, and 2'iTqR gives the connection length along the magnetic field lines. 
The parameters governing core or edge turbulence result from competition between ExB turbulence and parallel 
dynamics, and so the scale ratios enter the ion inertia and curvature parameters according to 

—J "^ = — 

respectively. These lead to the parameters governing the parallel electron dynamics, 

13 = Pee M=Tr^ g = M (39) 

These are the drift Alfven parameter, the electron inertia parameter, and the drift wave coUisionality parameter, 
respectively. For core turbulence C and jl are small, and (3 can be in the range from small values to several times 
0.1 in high performance tokamaks. For edge turbulence C and (i are larger than unity, and the turbulence becomes 
electromagnetic for (3 near or larger than unity. Ideal and resistive ballooning regimes occur when Pojb > 1 or 
Cojb > 1, respectively. The magnetic shear parameter is s, nominally given by dlogq/dlogr where r is the minor 
radius but generalisable to arbitrary geometry (cf. Ref. ^3)- usually of order unity. Under the shifted metric 

fluxtube geometry, s enters only through the shifts incurred in the y-coordinate while taking derivatives with respect 
to s ■ Further details on the significance of these parameters may be found in Ref. Q . 



III. ENERGETIC PROBLEMS WITH THE STANDARD GYROFLUID MODEL 



Most of the problems with the standard gyrofluid model 0, arise from the finite Larmor radius (FLR) terms. 
The closure approximations are done term by term in a reasonable way, but with subtle differences in the polarisation 
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equation and in the gyrofluid moment variable equations. In the polarisation equation the gyroaveraging of the 
distribution function is computed from density and temperature moments of a perturbed Maxwellian approximation. 
In the moment equations the gyroaveraging is done on the potential (concentrating here upon 0), which is subject to 
derivatives and then to the moment integrals. These two procedures are in general different unless one approaches 
them simultaneously with energy conservation in mind. For this reason, the result that the model has inconsistencies 
is not so unreasonable, particularly considering that its motivation starts with linear theory and only secondarily adds 
the nonlinearities one needs for turbulence. Additional to this is an inconsistency in treating the higher moments 
which arise from the curvature terms: a perturbed Maxwellian model is used for those terms while the model itself 
retains the parallel heat flux moments as dynamical variables. These difficulties are treated in turn, with the FLR 
terms first. The standard gyrofluid model is sufficiently well constructed that these repairs are in the end a minor 
matter. 

There are three gyroaveraging operators used in the model, given by Eqs. (20,26,27) of Ref. Q, 

4>G^Tn \^'J>o^h^^ %^^^6^(5ry^)^ (40) 

assuming a field aligned coordinate system and a Fourier representation in the two perpendicular coordinates, such 
that h = k\p1 (the model concentrates on a single ion species and leaves the electrons to be adiabatic). We relabel 
these in terms of Fi, F2, and F3, given by 

r.^rr r..*§ r3.^^,i.r.) (4i, 

and recast the gyroaveraged and FLR corrected potentials as 

4>G = Ti^ f2G = r20 f^G = r30 (42) 

for clarity. Note the factor of two inserted into the definition for Fs, so that VIq and VIq both have the same low-/c^ 
limit. 

The polarisation equation |17| is given by Eqs. (7,48), rewritten as Eq. (93), all from Ref. [J, 



1 + 5/2 (1 + 5/2)2 - 

preserving the normalisation of in terms of Te, where Ti_L is the perpendicular temperature disturbance defined as 
the normalised {v\ — 1) moment of the perturbed distribution function. Hence three gyroaveraging operators appear 
in the moment equations, but only two appear in the polarisation. Moreover, the Fade approximate forms are used 
in the latter but not in the former. Even if this were repaired, using 

in the gyroaveraging of the potential (cf. Sec. III.C.4 of Ref. 3]), it would be impossible to reconcile the fact that the 
operator does not appear in the polarisation equation, which in terms of the gyroaveraging operators reads 

ne = Fin, + F2r,± + ^^^^ (45) 
n 

For the purposes of energetic consistency it does not matter how or whether the Fg — 1 screening term is approximated, 
only which operators appear in the Ui and Ti terms. 
The ExB energy is given by 

^A^f (46) 
Using Eq. (|45l) we may replace this according to 

dA('-^l + ^]^ fdA^(^-^^l±^] (47) 
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Using the Hermitian property of the gyroaveraging operators and the definitions of the potentials, and placing the 
electron contribution on the right side, we find 

dA^_I^^^ [aA(l£h±Mi^^i::i\ (48) 



The drift energy is therefore expressed as a combination of potentials multiplied by thermal state variables. 

If the electrons are adiabatic due to fast parallel dynamics on closed flux surfaces, the electron density is itself 
replaced by the potential according to 

(49) 

where the angle brackets denote the flux surface ( "zonal" ) average. In this case the electron contribution is a sort of 
field energy which is combined with the proper drift energy to form the potential energy given by 




(50) 

The flux surface average is subtracted because the adiabatic state arises through the large value of the parallel 
wavenumber fc|| combined with the electron thermal velocity, in comparison with the dynamical frequencies, while for 
the zonal component there is no action by V||. The adiabatic electron approximation does not hold in general, but it 
is assumed in Refs. [^Q, and keeping it within this Section makes this discussion more transparent. 

The thermal free energy is given by the fluid moment variables in quadratic combination. The state variable free 
energy is given by 

r / n2 + T,-2 + f.2 \ 

(51) 




while the flux variable free energy (the "generalised parallel kinetic energy" ) is given by 

(~2 I ~ 2 I 2 ~ 2 \ 

Comparison of all these pieces leads to the conclusion that the combinations r^n^ + 0g and TiTi± + should appear 
together under first derivatives in linear terms in the moment equations in order that in combination among all the 
energy pieces the various terms reduce to terms involving single first derivative operators acting upon combinations 
of the variables, which have the form of total divergences, either SV|| or /C. It is clear that since no operation by the 
third gyroaverage T3 appears in the polarisation equation with only densities and temperatures kept in the closure 

for the total space density, there is no place for the third gyroreduced potential flc in the moment equations. This 
is the first and most obvious inconsistency of the standard gyrofluid model, and it impacts the parallel dynamics, 
the curvature terms (quasistatic compressible part of the drift dynamics), and the magnetic pumping process, each 
of which we presently examine in turn. Finally, the additional inconsistency in the treatment of higher moments in 
the curvature terms in the equations for g^n and qi± is addressed. 

We first look at the parallel dynamics. This involves conservative energy transfer in the sound waves and conductive 
heat fluxes. The simplest case of a sound wave in the absence of any effects due to the potential is of a parallel gradient 
in the pressure causing a parallel flow, and the corresponding parallel compression of that flow acting to restore the 
pressure disturbance. The total divergence term expressing energy conservation in such a local model as this one is 
i3V|| (pu||/_B), up to numerical constant factors dependent on how the temperatures are described. This divergence 
represents a transport process, in this case advection of thermal energy by u\\ . When the potential is also involved, 
charge currents become part of this, but the structure is the same. 

In the standard gyrofluid model, the part of the dynamics involving sound waves is 

drii uh ^ ^ 



dun 

— - = -Vii 
dt II 



(54) 
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in addition to the polarisation equation. Evolution of the potential energy, the thermal state variable energy, and the 
thermal flux variable energy pieces involved in this is given by 



'~dt 



^SV[|-^-^ + ii||V||0G 



dun 



«ll^ = -"itV|| 



tt,T, 



at 



-BV|| 



B 

n yrii + T,;i 

TiHiUn 



(^Hi + T,;||^ 



B 
B 



- -I- U||V||Tini 
+ ?I||V||TiTi|| 



(56) 
(57) 
(58) 
(59) 



When the pieces are summed, the transfer terms denoted by U|| V|| cancel, leaving the total transport divergence term 
denoted by 



-BVn 









Ti (rti + Ti [ 


) +0G 


^11 



B 



(60) 



that is, just the divergence of a transport flux given by the force potential in the equation for u\\ multiplied by 

W||B/i3. Here we note that Ti±^ is not involved, so there is no action by We find by this analysis that the sound 
wave dynamics is in order and does not require modification. The same conclusion results from examination of the 
non-closure part of the heat conduction dynamics, for which the equation parts are 



2 dt 



B 



and 



9% it 3 



TiTi\ 



'Ji± 



dt 



at 

-V|| 



-BV 



qi±_ 

B 



TiTi± 



and the energetics parts are 



1 

-TiTi 



dt 



-B\I\\ 
2. 



TiTi\\qi 
B~ 



and 



T,Ti 



j_L- 



dt 



= -BV 



3*^' II dt 



Vic 



+ '7j||V|| 
-9^11^11 



TiTi\\ 

TiTiw 



B 

qi± 



dt 



- + gi±V|| 

= -gij_V|| 



TiTi±_ - 
TiTi±_ 



no 



(61) 
(62) 

(63) 
(64) 

(65) 
(66) 

(67) 
(68) 



Here we note the factor of two difference in the definition of g^n here (conformal with the Braginskii definition ^^^) 
and in the standard model. 

There are also magnetic pumping terms in the gyrofiuid parallel dynamics, due to the combination of parallel fiow 
and conduction, magnetic moment conservation at the gyrokinetic level, and the parallel gradient in the strength of 
the magnetic field. Here, the standard model has the terms in the right places except for a single occurrence of the 

"forbidden" potential flc- 



du\\ 
'dt 



(t,x-T,||) 



fir 



V|| \ogB 



(69) 
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1 5Ti|| ^ . 



dqi 



dt 



dt 



= V||logB (71) 

VnlogS (72) 



If we merely replace flc by fie we restore consistency, 

dqi± 



n (Th, - Till) + Qg] V|| logs (73) 



dt 

When the free energy pieces are constructed, these effects form transfer channels which then properly conserve energy 

because TiTi± and CIq occur in combination. The reason is "forbidden" is that it doesn't appear in the polarisation 
equation. If it is present, then the conservation in the exchange between potential and thermal energy is broken. 
A similar problem appears in the curvature terms. For the thermal state variables in the standard model, these are 

dfj^ ^ ^ ^ + Og + TiPi^ ^ STjfj^ + Og + Z^G I ^^g^ 

Again, we need merely replace flc by Cla to restore consistency, 

^=^1^ 2 +^ 2 ) (^^^ 

so that TiTi± and fig again occur in combination. Once more, thermal free energy was already conserved in the 

standard model, but due to Slg a mismatch in the FLR part of the transfer between potential energy and thermal 
free energy remained. 

The curvature terms in the fluid moment flux variables present a different problem, the only inconsistency in the 
standard model which is not a FLR effect. There are no curvature terms involving the potential in the equations for 
U||, 5i||. find qi±. but there is a closure treatment at the level of the fifth moments which appears in the equations for 
the third moments (qi|| and qi±), arising from the factors of and in the grad-B and curvature drift terms in 
the gyrokinetic equation. In the standard model the closure for the 4th and 5th moments is taken from a perturbed 
Maxwellian (cf. its Eqs. 81 and 82). However, the model retains qi|| and qij_ as dynamical variables, and so the 5th 
moment should include contributions from pressures times conductive heat fluxes. For example, parts of the v^v^^ 
moment is provided under drift ordering by Pi±qi^\ and Pi]\qi±, and part of the f^f || moment is provided by Pi±qi±. 
In the normalisation, the factors of pi|| and pi± are replaced by unity. A way to do this systematically is to express 
the perturbed distribution function as a general six-term polynomial in which each of the coefficients is represented 
by one of the fluid moment variables retained in the six-moment model. Then, the fifth moments are computed by 
evaluating the integrals over w^W|| or ujj times the perturbed distribution function. The result of this calculation is a 
combination which automatically conserves free energy within the fluid moment system: 



du 



dt =§^(4^f||+25i||+5i±) (78) 
^ = |/C(3Sf||+8%||) (79) 
^ = ^/C(«||+6g,^) (80) 
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The "diagonal" terms in the imphed curvature matrix conserve automatically, but what this procedure has done is to 
ensure that the "cross" terms also act conservatively. 

The final consideration is the ExB advection terms, especially their FLR generalisations. All three potentials {(f>G, 

r^G, and flc) are involved in the standard model, acting through their respective drift velocities, 

ue^-F- V(^g we^-F- Vf2G We = -F • X/ha (81) 

where the tensor operator — F • V represents the more familiar (c/i3^)BxV. The standard model has "diagonal" 
terms, in which each variable is acted upon by u^; • V in its own equation, "cross" terms in which the pairs (n,;,Ti^) 
and are coupled by w^; • V, and then extra "diagonal" FLR effects in which Ti±_ and qi± are acted upon 

by 2W£; • V in their own equations. Because ilc is not ftc, the FLR potential energy VtGTi±/2 is not conserved by 

W^; • VTij^. Again, this is repaired by simply replacing Q,q by Vtc- 

Two minor considerations remain: first, the part of the collision process which enforces isotropisation should have 
the collision frequency multiplying Tn_-\-VLQ — Ti\\ rather than simply i so that the total dissipation is positive 

definite. This merely captures the correspondence between a gyrofiuid Ti^-\-VLQ and a fiuid Tj_L as discussed by Belova 
[l^ . Second, we have discussed energetics independently of nondissipative closure. The standard model constructs 
a curvature dissipation matrix in order to capture toroidal drift phase mixing by the velocity-dependent grad-B and 
curvature drifts. Part of this is nondissipative and can in principle capture the 5th moment effects discussed above. 
However, we choose here to separate these effects because in some applications involving nonpcriodic drifts there 
should be no phase mixing effect at all. Moreover, the dissipation matrix treatment itself is not a real success. 
It is explained (pp. 4057-8 of Ref. "^I) that a different set of coefficients was required to accurately represent the 
kinetic result close to marginal stability for the adiabatic-electron, toroidal ITG mode. This is in itself an admission 
of failure for the project of using the gyrofiuid system as a quantitatively exact representation of the gyrokinetic 
one. We do not attempt such an ambitious goal here; rather, the GEM model is intended for qualitative study of 
basic physics mechanisms, especially energy transfer between small scale turbulence and large scale fiows and MHD 
processes. Hence the neglect of dissipation free closure effects, and in general the neglect of dissipative effects other 
than collisions and Landau damping. 

The main point of this Section has been to highlight the way in which potential energy and thermal free energy are 
conservatively exchanged by the various mechanisms involved in low frequency ffuid drift dynamics and their capture 
by a gyroffuid model which is at least well behaved to arbitrary order in the finite Larmor radius parameter k±pi. If a 
turbulence model is to act at arbitrary FLR order it should conserve energy properly, even if for no other reason than 
that a numerical computation should not experience trouble in the spectral region around k±pi 1. The procedure 
to repair the standard model accordingly is to restore the use of the same closure treatment (and to the same depth in 
the number of moments kept) in the polarisation equation(s) and in the fluid moment equations. In the end, the same 
gyroaveraging operators appear in the polarisation and fluid moment equations — since the polarisation equation in 
the six-moment model only involves Ui and Ti± it only involves the flrst two operators Fi and and hence only two 
potentials (one nominal, c/jq, and one FLR, flc) can appear in the fiuid moment equations. The third gyroaverage 
operator may have a reasonable role in a gyrofiuid model which retains 4th and 5th moments (one level in the 
hierarchy past temperatures and conductive heat fiuxes), but we do not pursue this extension herein. 



IV. CONSTRUCTION OF THE GEM MODEL 



Subject to the conventions in Section^ our starting point is the polarisation equation, which links the variables 
riz and Tz± for each species to the electrostatic potential, 0. We neglect true space charge effects, setting the Debye 
length to zero and assuming the space charge densities all add up to zero. This is the statement of quasineutrality: 



Fn- 1 



2J-z± 







where the gyroaveraging and screening operators are defined separately for each species, 

To = Fo(fe,) 



(82) 



(83) 



with argument fe^ = k'^pl and squared gyroradius — /i^r^/B^, where B is the normalised strength of the equilibrium 
magnetic field. The considerations which lead to these forms are the ones given in the standard gyrofiuid model 0, 
as outlined in Section UTTl 
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We identify the generalised ExB energy using the polarisation densities in Eq. as 

^ To - 1 ^2 



Ue 



(84) 



Using Eq. 1821 and the Hermitian property of the F operators, we recast this as 



Up 



(85) 



where the gyroreduced potentials are given by 

= Ti^ 



To 



(86) 



Note that these are defined separately for each species and that there is no V^. 

We identify the thermal state variable part of the energy the same way as in the fluid models, 



(1/2)T. 



II ^ -'-^-L 



(87) 



The flux variable part of the energy is 



(2/3)g4 



9zi 



(88) 



We note here that since the evolution of Ut and Ue ultimately follows from the same moment equations, the combi- 
nations which must appear together under the V|| and K. operators in those equations are (pQ A-TzUz and -\-TzTz±. 
Observing this will guarantee energetic consistency. 

The magnetic field disturbances arise from the parallel magnetic potential, which is given by Ampere's law in terms 
of the total electric current, 



(89) 



The magnetic energy is given by 



which using Eq. may be rewritten as 



Ur, 



A 



^11 -^11 



(90) 



(91) 



This model neglects gyroaveraging and gyroscreening of the magnetic potential. On the same footing as the potential 
equation we would have a more general Ampere's law in which Uz\\ might be replaced by riM^n +T2qz±, and due to the 

same consistency considerations as with (jjc and flc several new finite gyroradius terms would appear in the magnetic 
flutter dynamics. This is being left for future work, however, because the consequences for energy conservation have 
not yet been worked out. 

In situations wherein the finite electron gyroradius is neglected, we simply have (f)G = 4' and Qq — for the 
electrons, so that the polarisation equation becomes 



Eirij + T2Ti± + 



Tn-l 



(92) 



where the species label is changed to i as it refers to the ions only. If as in most practical applications one takes a 
single component plasma with singly charged ions, we merely have Ui — 1 and r,; — Ti/T^, and with normalisation to 
that particular ion's mass, fii = 1 and /ig = —m^/Mi along with = t^, — —1. Since these are trivial restrictions, we 
present the GEM model in terms of the general forms following from Eq. H82|) . 
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The ExB advection operators and the nonhnear part of the parallel gradient are given in terms of Poisson bracket 
structures, 

U£-V=[0G,] we-V^Pg,] b^.V = -/3e[I||,] (93) 

in the two perpendicular coordinates, that is, 

^■^'^^ \dxdy dxdyj ^ ' 

Defined in this fashion, the quantities u^;, w^;, and bj^ are all divergence free as written; the generally finite divergences 
of u^; and are treated separately, through /C((/)g) and JC^fla), respectively. The advective time derivative and the 
parallel derivative are given by 

|^|+u..V V„=lB.V + b,.V (95) 

where B and B are defined in terms of the equilibrium magnetic field. The variation of B along a field line (poloidally 
around the fiux surface in a tokamak) is incorporated into the gyroradii 

Pz = (96) 
The perpendicular parts of the Laplacian are given by 

where g is the determinant of the {g^iv} elements of the entire metric, and is the perpendicular metric which 
involves only the two perpendicular coordinates {x and y). In these terms, the argument bz appearing in the T 
operators is given by 

hz = pI iKaTku) (98) 

where /i and v are summed over the two perpendicular coordinates only. Note that commutes with and that 
no gyroradius appears with the V^^ operator in the Ampere's law fEa. I89|l . 
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The six moment equations are the same ones appearing in the toroidal version of the standard gyrofluid model |j] , 
corrected according to the findings in Section Hill With and (fic, and Tz± and ^Iq, appearing together, and the 
background gradient forcing terms Uz and Tz (functions of x only) displayed explicitly, the dissipation free part of the 
equations for each species are 



dt 



We ■ V 



Tz + Tz 



Uzl\ 



TzPzll + TzPz± + 



(99) 



9^11 dUz\\ _ /- r _ n 

Pe-^ + IJ-z—JZ- + Pz^E ■ ^qz± = -V|| ( + [Pz + Pz\\\ 



dt 



dt 



4:Uz\\ + 2g^[| +qz± 



nG+TzTz±-TzTz\\]\/ll logs 



(100) 



dt 



Uz\\ + qz\\ 



= IC + TzTz« - (uzn + qz±) V|| logB 



(101) 



d Tz+Tzi_ 



dt 



Wfi • V ([ti^ + nz] + 2 



Tz + T, 



z± 



1)0 + ^0 + rzPz± + TzTz± 



BV 



qz±_ 
B 



gz±) Vii logs 



(102) 



dqz 



dt 



Tz+Tz\ 



IC flzTz 



3u 



Z\\ +092 II 



(103) 



t^z- 



dqz 



dt 



fiz^E ■ V {uz\\ + 2qzj_) = -V|| (flc + Tz 



Tz + Tzi_ 



IC [ M.t /"" '^f'^"^ ]~Tz(nG+ TzTzl_ ~ TzTzW ) V|| logB 



(104) 



The pressures are defined as 

Pz\\=nz+Tz\\ Pz±=nz+Tz± pz^Uz + Tz (105) 

and all the thermal state variables are operated upon by u^; • V, w^; • V, and Vii together with their gradients. 
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A. Local and Global Models 



The GEM model is variously cast in both local and global versions. In the local version the gradient drive terms 
appear explicitly in the equations as displayed above. Following the ordering, the pressures add linearly also in 
these quantities. The densities and temperatures are given prescribed forms, usually simple linear gradients, e.g., 
Uz = —i^nX such that w„ = |Lj_Vlogn2| gives the inverse of the normalised scale length, but they can be given 
arbitrary form, allowing computations within any prescribed gradients. For this version all dependent variables are 
given Dirichlet boundary conditions in the x-direction, 

/ = at ±^ (106) 

where is the domain length. 

In the global version the separate gradient terms do not appear in the moment equations. Instead, the fluctuat- 
ing gradient is part of the dependent variable. Accordingly, the dependent variables are given Neumann/Dirichlet 
boundary conditions in the a;-direction, 

^=0 at x = -^ / = at x^^ (107) 

dx 2 ■' 2 ^ ' 

allowing arbitrary profile evolution. In this version the curvature operator and also all the dissipation operators act 
upon the entire variable, including the profile. Consequently, the two dimensional equilibrium including parallel flows 
and currents, and heat fluxes, is also solved for and evolved self consistently with the turbulence. The transport 
problem is also solved. One usually operates in one of two limits: the run time is either longer or much shorter than 
the confinement time. In the latter case sources are not necessary; the profile is allowed to relax but is expected 
not to do so very much (perhaps 30% relaxation is acceptable). For edge turbulence in a thin radial layer such that 
Lx < L±, the confinement time is usually shorter than the time for the zonal flows, resulting from the evolving flux 
surface ( "zonal" ) averaged potential, to reach statistical equilibrium. In this case sources are necessary for all the 
state variables (n^ and Tz\\ and Tz± for each species). In either global case, the proflles cannot be prescribed except 
for the fact that the system is initialised with the proflles set into the state variables. 

A flnal consideration in the global model is polarisation: vorticity is given in general by the "gyrocenter charge 
density" made up by X]^'^^^^' noting that it is the total charge density that is set to zero. Here, the profiles are 
included in the dependent variables, so especially in the adiabatic electron version one must keep the unchanging (n^) 
in polarisation, 

(108) 

where the angle brackets denote the zonal average. In the local version, the profile function ni{x) does not appear in 
polarisation, since it is expected to be equal to ne{x). 

Both global and local versions are set up with globally consistent boundary conditions parallel to the background 
magnetic field. This ensures individual fulfillment of the periodicity constraint for each Fourier component in y as if on 
the entire flux surface, even if the toroidal mode spectrum is truncated lisj. Here we note that ky follows the toroidal 
mode number generally and hence the y-domain is periodic. The domain length for y is Ly. The periodicity constraint 
applies as a boundary condition in s after one single poloidal cycle, so the s-domain is always one connection length 
— TT < s < TT. Additionally, the j/-coordinate is shifted on each drift plane (constant-s surface), so that perpendicular 
dynamics is always computed with an orthogonal metric, and magnetic shear enters as a set of relative shifts in 
the y-coordinate in the expressions for d/ds, as explained in Ref. '20'|. This combination is required for capture of 
slab-character modes, of which the most important in turbulence is the nonlinear drift wave instability |7[l9ll2l|. 
Global consistency is also required to obtain the correct spectrum of sideband modes in the equilibrium, to which the 
turbulence and zonal flows are coupled by toroidal compression of the ExB velocity |22j . 



B. Fourier and Pade Versions 



The operators involving V^^, including the F's, are solved variously in xy-space or in kj^-space. The Fourier 
versions are set up to be compatible with either the local or global boundary conditions. The local model with 
Dirichlet boundaries uses half-wave Fourier transforms, with the basic x-domain odd-reflected about x — that 
is, f{Lx/2 + x') — ~f{Lx/2 — x') for < a;' < L^- The Fourier transform and its inverse is then applied to the 
doubled domain. In a similar manner, the global model uses quarter-wave Fourier transforms, with four copies of 
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the a;-domain arrayed odd-even-even-odd, and with the transforms apphed to the quadrupled domain. This Fourier 
version of either the local or global model is used whenever the argument = k\p\ of any species is expected to 
take large values. The hallmark example of this is ETG turbulence (Sec. EJ, which involves the entire scale range 
between pi and pe in a single component plasma. If the ion moment variables are kept, the Fourier version should be 
used to obtain an accurate ion response, which becomes more and more "adiabatic" (in the sense that — > —TiUe) 
with increasing hi. In current implementations, the Fourier version always uses the full FLR form for every species, 
including electrons, regardless of the expected values of b^- 

For standard ITG or edge turbulence cases, where the scale range reaches down to but not below pi, the Fade 
version may be used. This approximates the F's by Q 

ro(M->(i-plvi)"' (109) 
ri(M-(i-^pM) ' (110) 

r2(&.) = ipM(i-^pM) 'ri(&.) (Ill) 

However, if the finite gyroradius effects of more than one species are taken into account, one must solve the combined 
screening operator given by '^z{azl'''z){XQ ~ 1) to find 4>. This is simple in k_L-space but complicated if using the Fade 
forms in configuration space, even if one of the coordinates is Fourier decomposed. With only two gyroradii to follow, 
however, one has the acceptable operation involving two successive Helmholtz solves. In current implementations, the 
Fade version is only used for single component plasma cases in which the electron FLR effects are neglected. The 
adiabatic ion model is only used with the Fade version. The importance of the Fade version is that it is the only one 
easily generalised to fully global geometry, wherein the metric coefficients depend on x. 



C. Dissipation in GEM — Collisionless 



The only true dissipation in the collisionless GEM model is phase mixing due to the kinetic resonances caused by 
the parallel transit dynamics, i.e., Landau damping. This is represented by direct dissipation upon the parallel heat 
flux variables, using a Landau damping operator defined by 

flL. = flLO (l - 0.125l^qi?Vf ) (112) 

with constant a^o nominally set to unity, where qR is the field line connection length divided by 27r jl5j . V = t^J pz 
is the normalised thermal speed of species z, and V| is generally the full nonlinear parallel Laplacian divergence 

operator i?V||(l/i?)V||. Both qR and are normalised in terms of L^. 

The GEM model does not employ a curvature drift dissipation model. The standard one did so 4J, but also noted 
in detail the problems it raised. For the reasons discussed in Section ITTl| therefore, it is chosen to omit this feature. 
Nonlinear FLR phase mixing Q is also left out of the standard model ^ , for similar reasons of tractability. 



D. Dissipation in GEM — Collisional 

The electromagnetic gyrofluid model finds a very useful application in tokamak edge turbulence and hence 
requires a treatment of collisions which will capture the collisional Braginskii fluid model jl^ in the appropriate limit 
of large collision frequency and short mean free path. The usual types of dissipation are resistivity and thermal 
conduction, which in a model treating both velocities and heat fluxes as dynamical variables amounts to applying a 
dissipation matrix to their combination. For an isotropic temperature the method used in the DALE Landau fluid 
model 8J is sufficient. The gyrofluid model, however, additionally includes temperature anisotropy. To combine these 
effects, a simple drift kinetic Chapman-Enskog procedure is used to find the dissipative corrections to a state with 
stationary state variables |2^, generalised to a bi-Maxwellian distribution (defined by T^n, and Tz± for each 
species) for the dissipation free part. A Lorentz collision operator is used, and then the Braginskii coefficients are 
substituted to capture the collisional limit. The resulting model is given by 

dAn duz\\ 

Pe^- + = \- Pel^e 



r]J\\ H ( 9e|l + qe± + aeJ\ 



(113) 
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2~dr 



G 



(114) 



dt 



(115) 



dqz\\ 
dt 



(5/2) 



(116) 



dq, 



zl^ 



dt 



(5/2), 



-IJ-zVz 



qz± - OAazJw ) - 1.28i^2 - 1.5qz±] 



(117) 



where and az are the thermal conduction and thermal force coefficients for each species, and is the resistivity 
coefficient. Herein, the thermal force is kept only for electrons {az = ae), while for ions it is zero. For a single 
component plasma with singly charged ions, the values of the coefficients are 



T] = 0.51 ae = 0.71 Ke = 3.2 Ki = 3.9 



(118) 



To recover an isotropic model we would add (l/2)Tz\\ + Tz± to form (3/2)rz and qz\\ + qzi_ to form qz\\. This yields 
the forms in Ref. 8]. Then, we could recover the Braginskii formula for q^w by neglecting the inertia! and Landau 
damping terms in its equation (i.e., neglect qz\\ except in the coUisional damping term), as explained in Ref. 0. 

The terms in the temperature equations and the ones with the factors of 1.28 represent relaxation of anisotropy, 
and the others represent resistive (77) and thermal conductive (k^) dissipation. Note the combination TzTz±_ + in 
the temperature equations; this is required to make the dissipation positive definite, following the same considerations 
as those concerning energy conservation resulting from the same combination under the Vii and IC operators. 



E. Nonlinear Dissipation and the Numerical Scheme in GEM 



One final dissipation mechanism remains to be considered, and in gradient driven turbulence it is often the most 
important one: nonlinear cascading to arbitrarily small scales [2ll |. This enters explicitly as an artificial diffusion term 
in each equation for computations using a dissipation free scheme to calculate the nonlinear advection terms. 

The energy cascade in drift wave turbulence is generally local in k 1 -space p3 |. and can proceed in either direc- 
tion following the properties of the various nonlinearities 0, |2^ l2a |^ When energy cascades to the scale of the 
computational grid (highest k± values), it must be removed somehow lest the spectra approach the unphysical forms 
representing the maximum entropy state p5l | of the discrete system. One must check to ensure that the grid dissipation 
rate is independent of the resolution, essentially the same statement contained in high Reynolds number turbulence 
(dissipation independent of the diffusion or viscosity coefficient). 

In the past, the predecessor of this model has used an upwind scheme (a slope limiting algorithm |27| integratiner all 
the dimensions together (28i |. from computational fluid dynamics) implemented as discussed elsewhere |l6j| . Herein, we 
employ an alternative finite difference scheme which also does not involve Fourier transforming or spectral operations 
and hence is applicable to situations forbidding such operations. The scheme has been used with success by Naulin on 
the fiuid drift Alfven model [29i. The first derivatives involved in Poisson bracket structures are evaluated with the 
second-order version of the Arakawa spatial discretisation |30| . The linear terms involving parallel dynamics (d/ds) 
and perpendicular compressibility (/C) are evaluated with standard second-order central differences. Direct dissipation 
terms (e.g., collision-based frictional damping of Jy or Landau-based damping of heat fiuxes) are evaluated directly. 
The entire right side is evaluated thereby once per time step, but using a third-order "stiffly stable" algorithm derived 
by Karniadakis et al, according to which the previous three time steps of the dependent variables and the right hand 
sides are used to get the new time step '3l] . Since the entire right hand side is used this way, an unsplit second-order 
accuracy is achieved. Finally, the artiflcial dissipation terms are applied separately, using the dependent variables at 
the now-previous timestep. The structure of the equations is given by 



dF 



(119) 



where F is the functional of the dependent variables / appearing under the (d/dt) operator in each equation, S is 
the right hand side of each equation, and D is the artificial dissipation operator in each equation. The structure of 
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the algorithm is given by 
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So = Sifo) (120) 

(121) 



Fi^Fi+ MD{fo) (122) 
(apply boundary conditions) (123) 
/i ^ (124) 

where the subscript '0' refers to the current timestep, '1' refers to the new timcstcp, and the negative ones refer to the 
previous timesteps, Ai is the timestep interval, "boundary conditions" refers to the loading of the guard cells at the 
computational boundary so that derivatives are computed normally during the evaluations of S and and the last 
step recovering /i from Fi refers to the solving of the polarisation equations to recover (p and A|| and the evaluation 

of the gyroaveraging operators Fi and r2 to get the gyroreduced potentials 4>g and VIq. For waves, this scheme is 
stable without the use of D, allaying the principal consideration which led to the upwind scheme in the first place 
[T6|. But for turbulence we require the use of D. 

It is important to note that the artificial dissipation must work in all three coordinates, not just the two perpendicu- 
lar ones. ExB advection is the main agent causing the direct cascade towards large wavenumbers in the gyrofluid state 
variables, mostly (for edge turbulence) but also Tz\\ and Tzi_ (almost solely, for core turbulence). It is important 
to note that this occurs not only in k^-space but also -space, simply due to the statistics js^l- The dissipation 
operators must therefore function for both /c^ and fc|. 

One might be tempted to apply D to the force potentials, e.g., riz — (^g instead of 7^2, but this has been found to 
damage the solution measurably. It is indeed important not to apply artificial dissipation directly to either or ^||, 
the main effect of that being to destroy the Alfven dynamics (for A:|) or medium to large scale vorticity (for /c^). This 
was the problem with the upwind scheme as the kinetic shear Alfven velocity is scale dependent the exact one 
could not be used in the flux splitting involved in the scheme, so the fastest one was used (otherwise, the scheme is 
unstable). For /3e < m^jMi the fastest wave (following va) is at the lowest and the smallest scales (highest /cj^) are 
dominated by coUisional dissipation anyway (since v,. > Cs/L±), so edge turbulence was not strongly impacted. For 
core turbulence, on the other hand, the fastest wave (following Vf.) is at the highest k± so that the large scale MHD 
response at the lowest k± is strongly dissipated with a sort of super-resistivity acting directly upon A|| rather than 
J|| . To avoid the same problem with schemes with explicitly applied artificial dissipation, it is important to avoid 

application of any of the artificial dissipation operators directly to 
In the xy-plane the operations are summarised by the statement 

Us • V ^ Us • V - i^_LV^ - !^I|V| (125) 

in each equation; that is, artificial dissipation in both the xy-plane and the s-direction is applied to whatever is 
advected by the gyroreduced ExB velocity. An alternative is a hyperdiffusion for the xy-plane, so that 

Us • V ub • V + Vii^^Vi - i^l|Vj| (126) 

is used. For models with variable B these should be respectively changed to 

ub • V Us • V - V • v±pI'SIi_ - V • (biy||b) • V (127) 

and 

ub • V ^ Uij • V + V\v^pISJ\ - V • (bi/||b) • V (128) 

with 

pI-^ (129) 

in normalised units, so that the property of positive definiteness is preserved. If variable resolution causes problems, 
then the metric elements in these forms should be replaced by their flux surface averages. 

Typical values of these dissipation coefficients are set depending on the physical situation; in general they must 
be set as small as possible. Full resolution is found when it can be shown the resulting grid dissipation rate (not 
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necessarily the answer for the transport fluxes) is independent of the dissipation parameters. A resolution study will 
generally not be done at a particular value of the coefficient; rather, the coefficient should be made smaller when the 
resolution is increased. A window of operation opens when it is subsequently found that the above criteria for full 
resolution is met. Tests on core turbulence with adiabatic electrons {ve — P — He = 0) find that v± as small as 10^^ 
is possible with resolutions of ft,^ = /ij, = 1 or 2 x pg. Edge turbulence {C,$,ji all unity or greater; cf. Section Hll and 
Ref. 0) requires hx = hy = 1 x or smaller to be able to reduce to as small as 3 x 10^^. With = hy = 2 x ps 
a value of = 0.1 can be required, and this is generally too large to allow the vorticity dynamics in the range 
0.5 < k±ps < 1 to function properly. This is due to the robust nonlinear action by v^; • Vng in edge turbulence 
|2lj |. Under these conditions the hyperdiffusion form is necessary to be able to reproduce the nonlinear drift wave 
instability. For cold ion models (t^ = 0) with no temperature dynamics, this instability can be reproduced with a 
resolution oi = hy = 2 and a hyperdiffusion of i>± = 0.01. 

The parallel dissipation coefficients are easier as they are only needed to contain the cascade in the parallel wavenum- 
ber fc|| by the nonlinear perpendicular dynamics. Values of i^n = 3 x 10~^ for both ions and electrons are found to 
be sufficient with hs = 27rqi?/16, and for the most important wavenumber range —2 < k\\qR < 2 these lead to small 
corrections to the physical dissipation rates. It has been found necessary to use the same coefficient for both ions and 
electrons, to avoid artificial charge separation which can have a large effect on the spectral region with k^^qR moderate 
and k±ps small. 

V. SELECTED COMPUTATIONAL RESULTS 

It is not the purpose of this paper to enter detailed study of any of the problems the GEM model is to be applied 
to; rather, the focus is upon the way energetics works in the model and to use that to assist consistent construction 
of the model. Nevertheless, it is useful to apply the model briefly herein to an elementary situation whose capture is 
important (kinetic shear Alfven wave damping JJ3, 34J), a well known set of computational results (the Cyclone ITG 
turbulence campaign ^), and a demonstration that electron driven turbulence at scales below the ion gyroradius can 
be addressed with a gyrofluid model ("ETG" H^)- The latter two cases will be treated in proper detail in the future. 
ETG turbulence has been treated with a fluid model before, but only with adiabatic ion models [s^js^jH^- Herein, 
we apply GEM directly and flnd the ETG dynamics occurring naturally at its native scales. 

A. Kinetic shear Alfven wave damping 

Shear Alfven waves are well known from MHD |33|, and the coUisionless kinetic counterpart (KALF) is also well 
known '4^ . It has already been shown the gyrokinetic model treats them properly [sj |^ . We now use the result to 
calibrate the model Landau damping coefficient for the electrons. 

We take a basic parameter case with P — 1 and p. — I and e — 18350 with both coUisionalities set to zero as a 
reference. The magnetic field is straight and homogeneous {g^^ = g^^ = B = 1 with s = and K. — 0) and there 
are no background gradients (ujn — (^t = '-^i — 0)- The ions are cold (r^ = hence pi =0). The Fade version of the 
local model is used. With the homogeneous situation, the profile functions are set to zero and the domain is periodic 
in both y and s. The perpendicular domain sizes are = Ly — 2t:/K, with K = 0.1. The parallel domain is one 
connection length. The initial state is the sinusoidal disturbance = 10~^ (1 + cos i^Tx) cos Ky cos s, with = rie- 
The grid was 32 x 32 x 16 in {x, y, s}. The values of $ and p, and the Landau damping model coefficient a^o, were 
varied between 0.1 and 10. Artificial dissipation (vj^ and i/||) was set to zero. 

The KALF dispersion relation, shown in Fig.Q] has the two standard asymptotic limits (3/ p^ 1 and fc^ <C 1 where 
for these cases k\ = (5/4)if^. This range is found with the nominal sweep in /3 for /3 ^ 1, and in the sweep in p 
for /i ^ 1, and is well captured by the GEM model as shown by the comparison to the kinetic result using the root 
finding method of Ref. 34]. When fi ~ p there is substantial thermal electron resonance. In this regime the GEM 
model shows a peak, but with the peak value and its location only approximately captured. In the sweep of a^o the 
damping rate was proportional to a^o only for a^o < 1- For clo > 1 the effect is to remove the parallel heat flux from 
the dynamics, which becomes more ideal; hence the damping rate falls again. The maximum is found for a^o = 1-6, 
which is close to the actual thermal resonance at VS- The damping rate varies within the interval 8.5 < — 10'^7l < 9.6 
for 1 < < 2. As a robust model in the absence of fltting for all possible cases, it appears to be sufficient to simply 
leave a^o = 1; and the model performs qualitatively well. 
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B. ITG turbulence 

The standard of core turbulence studies with adiabatic electrons is the Cyclone project, which benchmarked a series 
of models and computations against a particular case of hot ion coUisionless turbulence and transport |3l • With the 
only free energy source being the ion temperature gradient, this is called ITG turbulence. With adiabatic electrons 
taken as a model {(3 — p, — C — 0), the parameter set is given by 

Wi3 = 0.290 e = 93.4 s = 0.78 

= = Ti = 1 LUn = 0.321 (130) 

The artificial dissipation coefficients are vj^ = 0.01, using hyperdiffusion, and i>\\ — 0.001. The boundary dissipation 
coefficient was 1.0. For the magnetic field the simple circular tokamak model with globally consistent boundary 
conditions and shifted metric coordinate system is used as detailed in Ref. [2(il |. The potential is initialised with a 
random disturbance bath in x and y 7] and a parallel envelope following the field lines from s — ,20] . and an RMS 
amplitude of lO^'*. The Fade version of the global model is used. The basic profile is given by 

P^'^''^ = Y ' TJ ^^^^^ 

and then both densities are initialised with uJnPo + (j) and the ion temperatures (Ti\\,T.i±) with LUiTiPQ. The adiabatic 
form of the polarisation equation with profiles is used, with {ue) = tOnPoi^): noting that tOt has no role for this 
problem. The perpendicular domain sizes are = Ly = SOtt, roughly commensurate with the global tokamak 
dimensions of a/ ps = 192 and 2TTr/qps = 350 (where r = a/2). The parallel domain is one connection length. The 
grid was 128 x 128 x 16 in {x,y, s}. 

The four cases with oji = {0.8, 1.0, 1.5, 2.2} were taken. Normalisation of the transport level is to i„, following Ref. 

0, so that the transport flux in units of the nominal Qi = ^(O.STin + Ti±)v'^'^, is recast in terms of a transport 

coefficient by taking Xi — Qi/i^n |VT|). Each run begins in a linear growth phase, overshoots to a transport level in 

the vicinity of Xi — 10; B^nd then saturates with the zonal flow dynamics (the part of the ExB flow arising from (^<PJ) 

reaching statistical equilibrium only well after t = 1000. Runs were taken to t — 4000. The transport is displayed 
statistically, with a sample taken at intervals of At = 10 in the phase 1000 < t < 4000. Slow relaxation of the 
temperature profile fills out the transport scaling curve. For each sample, the flux and gradient were averaged over 
all grid nodes in the part of the domain with < a; < Lx/4: before their ratio was computed. 

The resulting transport curve is shown in Fig.|21 wherein the triangle markers denote each sample (1200 in all), and 
the dashed curve is the fit to the gyrokinetic particle model results as given in Ref. H. Agreement at the 20% level 
is found for most of the curve, and moreover the nonlinear threshold agrees within the statistical scatter. Moreover, 
the fact that the groups of points from four decaying runs overlap well indicates the transport to be temporally local. 

C. ETG turbulence 

A class of turbulent dynamics at scales smaller than pi driven by VTg is called ETG |35|. Neither Vn nor VT!; 
is available as a drive because the ions are adiabatic (in the simplest treatments) or nearly so. Here, we carry both 
electrons and ions with the full six moments and allow the spatial scale range kept in the particular case to determine 
the dynamics. The Fourier version of the local model is used, with profile functions rig = = —uJnX and Tg = —ujtx 
and Ti = —ujiX. The same magnetic geometry as in the ITG examples above is used. The same random bath as above 
is used initially, but for rie- 

Here we merely demonstrate the ability of the GEM model to capture this ETG turbulence for typical core pa- 
rameters, the same as the one used for the ITG examples above, additionally with (3 — 0.464 and jl — 0.0254 and 
Ve = 0.0333. The artificial dissipation coefficients were = 3 x 10~^ (using simple diffusion, not hyperdiffusion) and 
v\\ — 10^"'. The initial RMS amplitude for was ao = 3 x 10"'^. The spatial domain size was = Ly = 47r/3 for 
the drift plane and one connection length along the magnetic field. The grid was 128 x 128 x 16 in {x,y,s}. The 
timestep was 5 x 10^*^. The run was carried for 2QL±/cs. With the minimum value of kyPi of 1.5, ITG activity is 
generally absent. The fastest growing spectral range is about 10 < ky < 20, representing structure scales Ay = n/ky 
on the order of lOpe- The spatial morphology shows these to be radially extended, with Ax > 8 Ay. The nonlinear 
transition begins at t « 8 and saturation occurs after t sa 12. A very strong transport level is found, just under 0.1, 
which in terms of electron scales is Xe ~ QptVe/ Lt- The spectrum is broader than in the linear phase, but still narrow 
compared to edge turbulence, and more importantly the transport spectrum peaks at ky = 10, very close to the linear 
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growth peak. The radially extended structures do persist in the saturated phase, by contrast to typical core ITG or 
edge turbulence. These features, shown for both linear and nonlinear phases in Fig. |21 are the same as those shown 
previously by nonlinear gyrokinetic studies • The amplitude of the electron moment variables rie and Tg y is about 
0.2, while the corresponding ion variables are about two orders of magnitude smaller. 

VI. SUMMARY 

The standard local gyrofluid model has been placed on energetically consistent grounds, with the moment variables 
and the electrostatic potential given a full finite Larmor radius (FLR) treatment at the same level of sophistication. 
The FLR effects on the magnetic potential and parallel velocities and heat fluxes is left to the future. With these 
changes it is possible to recover important results emerging from gyrokinetic computations, with a computationally 
more tractable model. Large systems may be treated with full resolution with what at present time may be regarded as 
modest computational resources. The model is flexible, to the extent that the level of sophistication can be increased 
or decreased while retaining energetic and geometric consistency. Both global and local situations can be treated. 
Highly detailed dissipative linear closures as discussed in the main references 0, 3 are not necessary in many cases 
of interest, in particular the one from the Cyclone study (Ref. 

For proper edge turbulence (/i > 1 and C > 1) in the electromagnetic regime (/3 > 1) the model functions much as 
in previous versions as published elsewhere Results from the two moment version GEMS (density and parallel 
velocity for both electrons and ions) are published elsewhere , showing the role of the three dimensional drift wave 
nonlinear instability in the context of tokamak edge turbulence as done previously pH for the corresponding fluid 
model. Work with cases with various ratios of rji = uJi/uOn (cf. l21*| for the role of this in the fluid model) and with 
two ion species is in progress. 
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APPENDIX A: SIMPLE CORRESPONDENCE BETWEEN FLUID AND GYROFLUID MODELS 



A simple exercise using the most basic reduced MHD interchange model helps gain insight into the relationship 
between the fluid and gyrofluid models. We start with the two equations written down in the Introduction, writing 
the electron density equation in terms of a charge density to equalise the units. For reasons which will become clear, 
we retain the diamagnetic compression effect in the density equation. We also incorporate the profile variation of 
the thermal state variables, normally acted up solely by ExB advection or magnetic flutter, into the corresponding 
dependent variables (the only difference this makes is that /C now acts upon the profiles, which is actually somewhat 
more realistic). The equations are 



B2 It 

dn, 



V\4> = -T,JC{n,) (Al) 



e-^ = n,e]C{4>) ~ T,JC{n,) (A2) 

with d/ dt representing the ExB advective derivative. We now define arbitrarily an auxiliary variable, N , as 

iVe = n,e-^^^Vi? (A3) 

The evolution equation for N is found therefore by subtracting Eqs. (Al,A2), 

dN ~ 

e— = n,eJC{(p) (A4) 
dt 

By inspection with any of the gyrofiuid models, we find that N is simply the gyrocenter ion density ti^, with the sole 
proviso that the background constant parameters for rie and rii are equal. This tells us that the MHD formulation for 
the ExB vorticity is identical to the cold-ion limit of the polarisation density in the gyrofluid model. The relation here 
is between the polarisation current in the fluid model and the polarisation density in the gyrofluid model. The MHD 
interchange term /C(ne) appears properly only if the diamagnetic compression effect is kept in the gyrofluid density 
equations (i.e., unhke for a fluid model, this effect cannot be neglected in a gyrofluid model). The ExB compression 
effect in the density cancels out of the interchange effect in the vorticity. However, it is necessary in either model to 
retain the ExB compression in order to conserve energy. The incidental benefit is to retain the geodesic curvature 
effect which is the principal mechanism limiting the growth of zonal flows p3 | . 

More detailed accounts of the correspondence between the fluid and gyrofluid models may be found elsewhere 

ET 
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FIGURE CAPTIONS 

Figure 1 . Kinetic shear Alfvcn damping rates versus normalised beta, electron mass, and Landau damping 
closure coefficient. In the leftmost two frames the blue line (whose peak is toward lower /3 values and higher /t values) 
gives the kinetic dispersion relation using the method of Ref. . The calibration works in the fig, 3> /ig regime and 
yields qualitatively similar behaviour elsewhere, though the details of the peaks can only be captured with a kinetic 
model. 

Figure 2. Transport of ITG turbulence found by GEM (triangles, one per sample as described in the text), 
compared to the gyrokinetic fit from Ref. 0. The transport diffusivity, Xii calculated temporally as described in the 
text, is normalised to a nominal value of xo — plcg/ Ln- The fact that the groups of points from four decaying runs 
overlap well indicates the transport to be temporally local. 

Figure 3 . Transport spectra and density morphology in core ETG turbulence in the linear (left) and saturated 
(right) phases, as described in the text. Lines marked 'n' and 'N' are for the particle flux, where it is positive or 
negative, respectively. Lines marked 't' and 'i' are for the electron and ion conductive heat fluxes, respectively. The 
scales are normalised to ps] multiply {x, y} and divide ky by 60.6 to obtain them in terms of pe- 
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FIG. 1: Kinetic shear Alfven damping rates versus normalised beta, electron mass, and Landau damping closure coefficient. 
In the leftmost two frames the blue line (whose peak is toward lower j3 values and higher jl values) gives the kinetic dispersion 
relation using the method of Ref. [3^ . The calibration works in the /3e ^ fie regime and yields qualitatively similar behaviour 
elsewhere, though the details of the peaks can only be captured with a kinetic model. 
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FIG. 2: Transport of ITG turbulence found by GEM (triangles, one per sample as described in the text), compared to the 
gyrokinetic fit from Ref. Ij. The transport diffusivity, Xiy calculated temporally as described in the text, is normalised to a 
nominal value of xo ~ p1ca/Ln. The fact that the groups of points from four decaying runs overlap well indicates the transport 
to be temporally local. 
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FIG. 3: Transport spectra and density morphology in core ETG turbulence in the linear (left) and saturated (right) phases, 
as described in the text. Lines marked 'n' and 'N' are for the particle flux, where it is positive or negative, respectively. Lines 
marked 't' and 'i' are for the electron and ion conductive heat fluxes, respectively. The scales are normalised to ps; multiply 
{x,y} and divide ky by 60.6 to obtain them in terms of pe. 



